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ABSTRACT 

We simulate the evolution of one-dimensional gravitating collisionless systems from non- 
equilibrium initial conditions, similar to the conditions that lead to the formation of dark- 
matter halos in three dimensions. As in the case of 3D halo formation we find that initially 
cold, nearly homogeneous particle distributions collapse to approach a final equilibrium state 
with a universal density profile. At small radii, this attractor exhibits a power-law behavior in 
density, p(x) oc |x|" rcrit , Tcrit - 0.47, slightly but significantly shallower than the value y - 1/2 
suggested previously. This state develops from the initial conditions through a process of 
phase mixing and violent relaxation. This process preserves the energy ranks of particles. 
By warming the initial conditions, we illustrate a cross-over from this power-law final state 
to a final state containing a homogeneous core. We further show that inhomogeneous but 
cold power-law initial conditions, with initial exponent y t > y cr i t , do not evolve toward the 
attractor but reach a final state that retains their original power-law behavior in the interior of 
the profile, indicating a bifurcation in the final state as a function of the initial exponent. Our 
results rely on a high-fidelity event-driven simulation technique. 
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1 INTRODUCTION 

The study of gravitating many-body systems is central to a wide va- 
riety of problems in astrophysics and cosmology. Perhaps the most 
important of these problems is the development of structure in the 
non-baryonic dark matter that dominates the mass budget of the 
universe. Dark matter is believed to be "cold," that is, the veloc- 
ity dispersion in the initial state is negligible. Such an initial state 
is gravitationally unstable, and subsequent evolution leads to com- 
plex three-dimensional structures such as halos, sheets, filaments 
and streams. Halos are the most recognizable structures. Structures 
such as filaments or streams inside halos are identified as coherent 
patterns in phase space, less readily visible but of significant im- 
portance. The phase-space structure of collapsed halos is important 
mainly because of its relevance for estimating rates and uncertain- 
ties for direct detection of particle dark matter (Maciejewski et al.| 
201 1 Vogelsberger & White 201 1 ) and for indirect detection from 



dark-matter annihilation near the Galactic center ([Hooper et al. 



2007 



2008 



Dobler et al. 2010 ) or in halo substructure ( Diemand et al. 



|The Fermi-LAT C ollabora tion: M. A ckermann et al.|2011) 



Currently, the dynamics of these systems can be studied only with 
large numerical Af-body simulations. 

Cosmological Af-body simulations reveal a number of com- 
mon features of halos that so far have resisted analytic explana- 
tion. The most striking feature is that the spherically averaged den- 
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sity profiles of collapsed halos appear to have the same form, apart 
from scaling factors for the total mass and radius. The common pro- 
file can be approximated by the Navarro, Frenk & White ( 1996 ) or 



Einasto & Einasto 



NFW formula, p(r) oc r _1 (r + a)~ 2 or thel] 
formula p(r) oc exp[-(r/a) a ]. A more detailed and better motivated 
analysis of shape universality is given by Dehnen & McLaughlin 
(2005|. See |Merritt et al.| ( [2006| ) for a comparison of models. The 
existence of a "universal" final state for equilibrium halos is rem- 
iniscent of Lynden-Bell's (1967) argument that chaotic changes in 
energy caused by the rapidly varying gravitational potential dur- 
ing collapse lead to an equilibrium state that can be determined 
by statistical mechanics. Lynden-Bell called this process "violent 
relaxation" and we shall use this term as well, although without 
the implication that the final equilibrium state can be computed by 
methods of statistical mechanics or that the changes in energy are 
chaotic. 

By their nature, Af-body simulations of fully three- 
dimensional systems probe only a limited range of length scales. 
Systems of lower dimensionality, or having special symmetries, 
provide laboratories in which to explore many-body dynamics with 
far higher phase-space resolution than is otherwise possible, and 
studies of these simpler systems can suggest explanations for be- 
havior seen in the three-dimensional case, such as the emergence 
of a universal halo profile (see |Dalal et al.| [2010 and references 
therein). Furthermore, simulations of reduced systems may ex- 
hibit some of the same systematic problems as simulations of 
three-dimensional dynamics, providing an opportunity for high- 
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fidelity tests of the overall methodology. Recently, Binney ( 2004 ) 
argued that simulations of one-dimensional systems indicate that 
discretization errors affect the central structure of simulated col- 
lapsed objects. He also observed that violent relaxation led to a 
power-law density profile, p(x) oc \x\~ 7 with 7 ~ 1/2, over three 
orders of magnitude in radius, although this was not the first time 
that this type of collapse had been investigated ( Gurevich & Zy- 
|bin|19 95). Although he examined only one set of initial conditions, 
Binney conjectured that this power law was universal over a wide 
variety of "cold" initial conditions. An understanding of this scal- 
ing could illuminate the origin of the central density cusps found in 
simulations of three-dimensional collapse. 

In this paper, we report simulations of the formation of cold 
dark matter halos having "slab" symmetry, in which the density 
is independent of two of the three spatial Cartesian coordinates. 
This system is equivalent to a one-dimensional system of gravi- 
tating particles, as studied by Binney and others. Our simulations 
improve the phase- space resolution of previous work by more than 
an order of magnitude, and explore a much wider range of initial 
conditions. We find strong evidence that cold collapse leads to final 
equilibria that are close to a universal power-law density profile, 
with an exponent less than 1/2. 

When the initial conditions are slightly "warmed", we ob- 
serve a cross-over from the power-law behavior of cold systems 
to the cored behavior displayed by warm systems. This dichotomy 
raises the question of the meaning of "coldness", which we dis- 
cuss and relate to our simulation results. We have also examined 
inhomogeneous initial states, with power-law density profiles, and 
observed an apparent bifurcation in the dynamical behavior of 
the system as a function of the initial power-law exponent. Ini- 
tial states with power-law density profiles shallower than |x| _7crit 
with y crit ^ 0.47 evolve to a universal profile |x| _7crit ; initial states 
with profiles steeper than |x| _7crit display a persistence of their ini- 
tial power law. Therefore, the |x| _rcrit profile represents a type of 
dynamical attractor, accessed by "shallow" initial states. "Steep" 
initial states are not attracted to the critical scaling solution. 

In section [2] we explain the basic setup for one-dimensional 
gravitational collapse and discuss the elementary properties of sta- 
tionary self-similar solutions. In section [3] we describe the algo- 
rithm used to perform the simulations in this work. Section[4]details 
our numerical results on the density profile, the particle energies, 
and the scaling properties of the final state. Section [5] summarizes 
our findings. 



2 ONE-DIMENSIONAL GRAVITATING N-BODY 
SYSTEMS 

One-dimensional gravitating (IDG) systems provide the simplest 
possible framework in which to investigate the nature and out- 
come of gravitational collapse. Violent relaxation and phase mix- 
ing ( [Binney & Trem aine 2008, Kandrup 1998} in such systems 
have long been a subject of computational study fCuperman et al.| 
[T9691 |Lecar & Cohenl[T97TT ). Simulations of IDG structure for- 
mation with cold dark matter and cosmological initial conditions 
date back to Doroshkevi ch et al.] ( |1980| >. Because of our interest in 
cosmological systems with very large numbers of particles, we fo- 
cus on the collisionless limit where the number of particles tends 
to infinity, holding total mass and energy fixed, and keeping time 
bounded ( Braun & Hepp 1977). The mean-field dynamics of an iso- 
lated system described by this limit satisfies the collisionless Boltz- 
mann equation (CBE), with potential and forces determined by the 



Poisson equation; 



d l +v d l 
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= 4nGp = AnGm p \ fdv; 



(1) 

(2) 



Here x and v are position and velocity coordinates in the two- 
dimensional one-particle phase space, f(x, v, t) is the particle num- 
ber density in this space, p{x, t) is the mass density, 0(x, t) is the 
gravitational potential, and m p is the particle mass. In this limit, 
two-body interactions are negligible. For collisional IDG systems 
see |Rybicki| ( fT97T}|Yawn & Miller| (1997} and |Gabrielli et ah] 
(2009) . 

If the potential is stationary, the specific energy E = ^u 2 +$>(x) 
is an integral of the motion. By Jeans' theorem, stationary solutions 
of the collisionless Boltzmann equation can only depend on the 
integrals of motion so f(x, v) = f(E). Then the Poisson equation 
(2} becomes 

M =SnGm »i (2£- • (3) 



We may multiply this by dO/dx and integrate once to obtain 

d ®\ 2 ™ i, dw » r r dE fW 

— = w(O) where — = 8/rG — — — -z. 
dx) dO J (E-<$>yi 2 

Re- writing this equation in the form 

dx _ 1 
dO ~ " Vvv(0)' 

we see that all stationary solutions must have left-right symmetry 
about some point x which we may choose to be the origin, i.e., 
O(x) and p(x) are even functions of x. 

Among all stationary solutions, a family of self- similar solu- 
tions can be constructed by choosing the phase- space distribution 
to be a power law in energy, 



(4) 



(5) 



f(E) = E- 1/2 (E/E r p . 



(6) 



We have chosen to make the position coordinate dimensionless, and 
factors of E have been inserted so that / has dimensions of inverse 
velocity. Computing the density from equation (5} gives 

p(x) = V5fl(I,p - I) m p (O/E ) l/2 - p , p>\, (7) 

where B{a, b) is the beta function. Choosing a center of symmetry 
where the potential is set to zero, we write the potential, density, 
and enclosed mass in the form 



p(x) = p \x\ r , m(x) - m \x\ 

where 

4p-2 



i-r 



OW = ® \x 



2-r . 



r = 



1+2/?' 



O(0), (8) 



(9) 



The constraint p > - implies 7 > 0, and the requirement that any 
compact interval containing the origin has finite mass implies 7 < 1 
or p < |. The relation between E , O , and p is straightforward to 
derive but is not needed here. 

Self- similar solutions can also be described in terms of the 
action, 



(10) 



\ r 2 r* JCmax 

J(E) = — (j)vdx = - J yj2E - 20(x) dx, 

where E = 0(x max ), and x max is the maximum excursion of the 
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orbit from the center. For a power-law stationary state, the relation 
between energy and action is 

2 3/2 F l/2 / F \ 1/(2-7) 

w -5iR>(£) B (*-»- 

Since / oc £(4-r)/(4-2y)^ me phase-space distribution function {6} has 
the power-law form 

f(J) OC /-( 2+ r)/(4-r) j-4p/V + 2p)^ (12) 

Because of its relation to energy, the maximum excursion of an 
orbit is also a useful surrogate for investigating the scaling of the 
potential. 

The emergence of a self- similar final state from cold collapse 
in one dimension was observed by Gurevich & Zybin (1995) and 
Binney ( 2004). An approximate analytic theory for the self- similar 
state was presented in |Gurevich & Z ybin (1995), with an estimated 
exponent y ^ 4/7; simulations described in that work were claimed 
to be consistent with this approximate exponent but few details 
were given. Their theoretical estimate is unlikely to be accurate 
because it is based on the assumption that adiabatic evolution be- 
gins at the transition from one to three streams; we have observed 
empirically that the actions evolve after this time. The exponent 
observed in Binney ( 2004) was estimated there to be y ^ 1 /2. 

The meaning of coldness in the initial states requires some dis- 
cussion. Intuitively, a state is cold if the velocity dispersion is small 
compared to the average potential. In the continuum limit, a per- 
fectly cold state corresponds to a phase- space density that is sup- 
ported on a smooth curve in the (x, v) plane, with only one velocity 
at any location. In the context of a particle simulation, a state is cold 
if it well-approximates such a continuum state. These statements 
fall short of being a precise definition, but the meaning of coldness 
will be further examined in the context of our simulation results, 
described below. A cold system in the continuum limit remains a 
sub-manifold in phase space indefinitely. However, a particle real- 
ization of a cold initial condition could depart from the behavior of 
the continuum model, either due to particles leaving the lower di- 
mensional surface or by under-resolving it. These effects can artifi- 
cially heat the system by introducing additional velocity dispersion 
to a coherant stream of particles at a given location. We will see 
that the low dimensionality of the phase- space surface is very well 
preserved in the particle realizations, until late times when discrete- 
ness effects grow to importance. 

Any coordinate parametrization for a curve in the (x, v) plane 
provides a Lagrangian coordinate for particles along the curve. Par- 
ticle number or rank along the curve is one such coordinate. Up 
to an irrelevant numerical factor, and allowing for continuation to 
negative values, the "mass enclosed on the curve" is equivalent to 
this particle number. In a self-similar state, this Lagrangian mass 
coordinate satisfies a power-law relation of the form 

fi(J) oc f f(f)df ocjQ-W^, (13) 
Jo 

using the expression from equation |T2] >. Note that the expression 
in equation (T2| is for equilibrium systems that are not necessarily 
cold, and that in this context is only valid at times when the system 
is very tightly wound and the orbital phases are mixed. Identifying 
jd with initial particle rank, this relation provides another avenue for 
investigation of the scaling properties of the final stationary state, as 
long as this identification remains valid. The identification becomes 
invalid when J is no longer a monotonic function of particle rank. 



3 ALGORITHM 

The system studied here consists of N parallel sheets, each of sur- 
face density H/N, interacting only through their mutual gravita- 
tional force; we shall call these sheets "particles" and treat them 
as an effective one-dimensional gravitating system. When two par- 
ticles collide, they can either be passed through each other with no 
change in velocity or scattered elastically; the resulting dynamics is 
the same in either case, up to a possible relabelling of particles. For 
definiteness, we choose to let particles pass through each other at 
collisions. With this choice, the initial positional rank of particles 
in a cold state remains a suitable Lagrangian coordinate throughout 
the evolution. 

The potential energy between two particles is 

<f>y =2nG(L/N) 2 \xi-Xj\ 9 (14) 

such that the force on particle i due to particle j is 

Fy = -d^tj/dxi = -2ttG (27A0 2 signte - xj). (15) 

The resulting net acceleration of particle i due to all other parti- 
cles is dependent only on its positional rank, k u with respect to the 
median position. Particles to the right of the median position have 
ki > and those to the left have k t < 0. Simulations with odd parti- 
cle number have a particle at x = and integer values of rank. Sim- 
ulations with an even number of particles have half-integer rank. 
Only particles j satisfying \kj\ < \k t \ contribute to the net force on 
the particle i; forces from \kj\ > \k t \ are balanced right and left. In 
the rest of the paper, and in all our reported simulations, we choose 
units such that 2nG = 1 and £ = 1 . 

The form of the force law implies a great simplification for 
simulations of this system. The acceleration remains constant as 
long as the particle's rank k { remains unchanged. Between particle 
crossings, the orbits Xt(t) are quadratic in time t, and the time of 
the next crossing of particles i and j can be obtained by solving the 
quadratic equation Xi(t) = Xj(t). Therefore, a numerical simulation 
can follow the exact evolution of the system, with errors arising 
only from roundoff in the solution of the quadratic equation and in 
the machine representation of phase- space coordinates. 

Conceptually, the calculation consists of the following steps. 
Compute all crossing times between each particle and its two neigh- 
bors on either side; find the earliest crossing time; advance to that 
time and update all particle positions; swap the crossing particles 
and adjust their accelerations; recur se. However, since the accelera- 
tion of non-crossing particles is unaffected by the crossing, it is not 
necessary to update the positions or re-calculate crossing times for 
non-crossing pairs or their neighbors. Nor is it necessary to advance 
non-crossing particles in time, except when the total particle state 
is synchronized for output of a data snapshot. Such simulations are 
event-driven rather than driven synchronously in time. The events 
of the simulation are particle collisions. Similar methods for com- 
puting the dynamics of one dimensional gravitating systems have 
been developed over the last two decades, greatly improving the 
reach of simulations jMineau et aT]|1990| |Yawn & MiTTer||1997| 
INoullez et aL|2003) . 

One possible code implementation uses a heap data structure, 
which stores particle pairs sorted according to their crossing time 
(Noul lez et al.||2003) . This reduces the cost of locating the next 
crossing to <9(log AT), which for very large N dominates the opera- 
tion count, since all other operations per crossing are 0(1). How- 
ever, with this approach it is not straightforward to account for 
changes in crossing times between particles i and i + 1 just crossed 
and their neighbours i - 1 and i + 2, since the heap position of these 
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adjacent pairs is unknown. We developed two different strategies, 
implemented in two completely independent codes, to overcome 
this problem and still preserve the 0(logN) operation count for 
finding the next crossing. The first method uses a search tree in- 
stead of a heap. A further speed-up is achieved by integrating posi- 
tional offsets between neighboring particles instead of the particles 
themselves. The second method allows the heap to contain dupli- 
cate collision sets, marked to indicate age of the sets. To suppress 
the accumulation of roundoff error, we found it important to use 
80-bit floating-point arithmetic (64-bit mantissa). Simulations run 
with ordinary double precision arithmetic (53-bit mantissa) show 
clear signs of roundoff errors, starting with a loss of local phase- 
coherence in simulations from cold initial conditions. 

All simulations presented in this paper explicitly enforce left- 
right symmetry about x = 0, implemented using a mirror bound- 
ary condition, with the particle state reflected about x = 0. This 
effectively doubles the resolution; the numbers N reported in the 
sections below reflect the total number, including mirror particles. 
By running simulations without an imposed symmetry, we verified 
that our conclusions do not depend on this enforced symmetry nor 
on the symmetry of the initial conditions. 



T= 1 8 T= 25 




-2 
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4 RESULTS 

4.1 Description of the Collapse 

The simplest cold initialization places particles on a smooth curve 
very near to but slightly displaced from the x-axis in the one- 
particle phase space. Placing particles with equal spacing in x gives 
a homogeneous initial density, and the velocities can be assigned 
as a smooth function of x. It can be shown analytically that initial 
conditions such as the Hubble flow with v oc x result in an endless 
repetition of homogeneous collapse and expansion, with no phase 
mixing or energy transfer among the particles. Initial velocities de- 
fined with a smoothly varying non-linear function of x, however, 
will generate a particle locus that forms an ever-tightening spiral 
curve in the (x, v) plane. This process is illustrated in Fig.[T] which 
shows four snapshots of the phase space evolution of our fiducial 
simulation, described below. We refer to the different arms of the 
phase space spiral at a given position as streams. The particles are 
divided into five quintiles in initial position, and these are color 
coded in Fig. [I] in a spectrum ranging from blue to red. Blue par- 
ticles began closest to the center |x| = 0, and red particles began 
farthest from the center. The density profile of the collapsed object, 
to be discussed later in detail, is the projection of the phase-space 
curve onto the x axis. Vertical portions of the phase-space curve in 
projection create caustic structures in the density profile (examples 
can be seen in Fig. [13]). 

Early in the evolution, rapid changes in the gravitational po- 
tential cause energy transfer among the particles, leading to violent 
relaxation (Lynden-Bell 1967). Note that the orbital specific ener- 
gies of particles, E k - \v 2 k + $>(x k , t), do not sum to the conserved 
total energy, E tot - \Y^k i v l + ®( x &>0]> so me energy changes of 
particles need not sum to zero. During this period, the potential os- 
cillates in time, with oscillations damping as the particles approach 
their final relaxed state. Violent relaxation ends when the particle 
locus is tightly wound and the gravitational potential is no longer 
changing significantly. 

Although the gravitational potential approaches a stationary 
state, the spiral curve describing the locus of the particles contin- 
ues to wind up. In the limit where the state is described by a con- 
tinuous curve in the (x, z/)-plane, this spiral will continue winding 



Figure 1. The evolution of an initially cold system of 10 4 particles with 
perturbed initial velocities. The simulation is described at the beginning of 
section |4~2] The spiral locus of the particles winds up as time elapses. The 
particle colors are arranged in a spectrum that encode the initial particle 
positions. Blue particles began the simulation closest to the center x = 0, 
while red ones began farthest from the center. The upper left panel shows 
two times at and near the start of the simulation. 

indefinitely. In the simulated system, with a finite number of par- 
ticles, the particle density along the spiral eventually decreases to 
the point where the curve is no longer identifiable; at this point, 
discreteness effects have erased any visual remnant of the phase 
coherence necessary to define such an approximate curve. 

These simple observations are directly related to the meaning 
of convergence in the collisionless limit. In order to make con- 
tact with the stationary solutions of the collisionless Boltzmann 
equation ([TJ, we must understand the nature of this convergence as 
N — > oo. A state supported on a smooth curve in phase space for- 
mally has an infinite number of particles, but it never becomes sta- 
tionary in a pointwise sense. Rather, in the collisionless limit, these 
states converge only weakly, as measures, to stationary solutions of 
the CBE (Braun & Hepp 1977). Therefore, the winding spiral locus 
can be identified with a stationary solution only through a process 
of coarse-graining. Coarse graining need not be restricted to aver- 
aging over phase space coordinates, but could also be implemented 
by averaging over orbital phase. For example, the / dependence in 
equation (T3] > relies on such a coarse graining. 

4.2 Cold initial conditions with homogeneous density 

4.2.1 The fiducial simulation 

To obtain the density profile for the final state of cold collapse, we 
ran a suite of simulations, with particle numbers clustering around 
N = 10 4 . The choice of N in these simulations is discussed be- 
low. For our fiducial simulation, we chose 2nG = 1 and set the 
total mass to one, £ = 1, implying a system of units where typical 
lengths and dynamical times are of order unity. The particles are 
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Figure 2. The density profile at a time t = 200 for a suite of cold simu- 
lations, with initial conditions as described at the start of section |4~2] The 
blue line is the density profile averaged over six simulations with slightly 
different particle numbers. The shaded region is the error on the mean as 
determined from the variance of the six simulations. The thin black line 
shows a power law x~ 1 ^ 2 . The red dashed line shows the shallower power 
law obtained by fitting the relation J(p) of equation fl3) . The thick black 
line shows the model in equation fTT) . 

initially equally spaced over the interval [-tt/2, tt/2], and the initial 
velocities are set by 

Vi(xj) = -Vo sin(xj), (16) 

with the scaling factor V = 0.001. This state represents an ini- 
tially overdense region near turnaround. The evolution of this state 
leads to a single collapsed object with a stationary profile. The sim- 
ulations were run to time t = 700, although several results are re- 
ported at time t = 200. At time 700 (200), the outermost particles 
have executed roughly 150 (40) orbits, with inner particles execut- 
ing orders of magnitude more. The simulation at these times has 
evaluated approximately 1.24 x 10 10 (3.5 x 10 9 ) particle crossings. 

Fig. [2] shows the density profile of the equilibrated system. 
The number density is computed using bins of variable width, set 
down so that there are approximately 20 particles per bin, and the 
abscissae for these points are placed at the geometric mean of the 
bin edges. The curve is the average of six different simulations. 
As is clear from Fig. [2] the final state is self- similar over four or- 
ders of magnitude in distance from the center. The fine- scale struc- 
ture in the average profile has two easily identifiable sources. First, 
the pointwise density profile for each simulation is dominated by 
caustics, the heights of which are regulated by the discreteness of 
the particle representation. Second, approximating the continuous 
density distribution along a curve in phase space by discrete parti- 
cles leads to fluctuations in density because of counting statistics. 
The caustic structure is easily identifiable at early times, but at the 
time represented in the figure, the discreteness of particle sampling 
makes it impossible to identify individual caustics; in all likelihood 
there are multiple caustics per bin. 

The density estimator for a single simulation approximates the 
continuum density profile but suffers from the fact that the bin spac- 
ing and the particle phases interact to produce slight (discretized) 
variations in the measured density in the bin. Roundoff error may 
also be responsible for occasionally pushing particles over the bin 
boundaries, a form of discretized phase error. To explore the impact 
of discretization and roundoff error, we have repeated the simula- 



tion with slightly different values of N; the six runs used to produce 
Fig.[2]have N = 10006, 10012, 10018, 10024, 10030 and 10036. 
The bin positions are the same for all six simulations, allowing us 
to quantify the sensitivity of the density estimator to the types of 
phase errors probed by these different choices. The mean density 
obtained from these six simulations is shown in Fig. [2] as a jagged 
blue line. The uncertainty in the mean, estimated from the variance 
of the six runs, is shown as a shaded region. These "errors" are 
generally somewhat smaller than Poisson in the power-law region. 

The fluctuations in the estimated density should not be inter- 
preted as Gaussian random fluctuations. An Anderson-Darling test 
run on the 6 values of the density in each bin suggests that only a 
small fraction of the bins, less than 15%, are consistent with hav- 
ing been drawn from a Gaussian distribution. A histogram of the 
skewness of the 6 density measurements in each bin shows no net 
skewness; however, most of the bins show significant negative kur- 
tosis (rounded peak with very short tails). Since fluctuations are 
principally generated by a small number of particles shifting across 
a single bin boundary due to phase differences, we would expect 
particle excursions to be at most one bin division and the tails of 
the distribution to be extremely small. We have checked that the 
distributions of skewness and kurtosis do not evolve as a function 
of distance; the inner decades display the same skewness and kur- 
tosis properties as the outer decades. Therefore, we feel confident 
that we can use the measured variances as weights when fitting a 
power law of the form p(x) = A|x| _r , without systematically over 
or under weighting any region of the data. However, it would be 
unwise to propagate errors from p(x) to errors on the fit parameters 
A and 7, because such errors would be very difficult to interpret. 
For this reason, we use the bootstrap technique to generate a large 
number of mock datasets, with replacement. We use a least- squares 
procedure to fit the parameters for each of the mock datasets, and 
report the mean and variance of these fits to the synthetic data. 

Using this method we have fit a power law p(x) oc \x\~ 7 to 
the portion of the data in the range 10~ 3 < \x\ < 10 _1 . We find 
y = 0.472 + 0.001 with a reduced x 2 ed - 0.339, computed using 
the standard deviations within each bin. The small value of the re- 
duced x 2 compared to the expected value x 2 ed - 1 arises from the 
negative kurtosis discussed earlier. The solid black line in Fig. |2]is 
obtained by multiplying this fit with a smooth truncation function 
(reminiscent of Einasto & Einasto 1972 models) 

PmodelO) a: \X\~ 7 exp (- (\X\/X ) 2 ~ 7 ) , (17) 

but the truncation of the power law in the outskirts is not the focus 
of our study. The fiducial simulation shows that the density follows 
a power law over at least four orders of magnitude in distance from 
the center, the extent of which is only limited by our finite resolu- 
tion. This confirms the results of Binney (2004 ) with many more 
particles (301 vs. 40000 in our highest N simulation). However, 
the slope y = 0.472 + 0.001 differs slightly but significantly from 
the value 1 /2 suggested by Binney, shown in Fig.|2]for comparison. 
There is a small systematic trend in the fit as a function of the range 
of data being fit; decreasing either the inner or the outer boundary 
causes y to be slightly smaller (shallower), though in every case 
significantly different from 1/2. 

4.2.2 Resolution study 

In Table[T] we show how these results depend on particle number N 
over the range 2 x 10 3 to 4 x 10 4 . Columns 2 and 3 show the best fit 
values of the normalization A and exponent y of the density profile 
at time t = 700, Column 4 shows the reduced ^ ed ; this should be 
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Figure 3. The best fit value of y, with errors from the bootstrap procedure described in section pT2T] The variability from one snapshot to the next is higher 
than the bootstrap errors would suggest. This is because the bin locations in the density estimator shift between snapshots, but do not shift in the bootstrap fit 
to a single snapshot. The time average value of y, and the rms fluctuation with time can be found in column 5 of Table[T] 



N A y x^d <7>t 7m binned y m single run y from J{ji) 

2000 0.3499 + 0.0302 0.4345 ± 0.0255 0.6307 0.4778 ± 0.0187 0.4699 ± 0.0009 0.4828 + 0.0003 0.4155 + 0.0008 

4000 0.3048 + 0.0065 0.4602 + 0.0043 0.2583 0.4710 + 0.0152 0.4602 + 0.0007 0.4499 ± 8e-05 0.4129 + 0.0004 

10000 0.3088 + 0.0043 0.4579 + 0.0031 0.3665 0.4703 + 0.0099 0.4567 + 0.0002 0.4378 ± le-05 0.4115 + 0.0002 

20000 0.2865 ± 0.0015 0.4720 ± 0.0014 0.3393 0.4722 ± 0.0089 0.4497 ± 0.0001 0.4478 ± 8e-06 0.41 10 ± 8e-05 

40000 0.2905 ± 0.0007 0.4700 ± 0.0006 0.3102 0.4691 ± 0.0068 0.4491 ± 7e-05 0.4606 ± 4e-06 0.4107 ± 4e-05 



Table 1. We examine the sensitivity of the fits to the number of particles being used in the simulation. Each line represents six simulations with slightly 
different particle number, with our fiducial set of initial conditions. Column 1 shows the approximate number of particles. Columns 2 and 3 are the best fit 
normalization and scaling for the functional form p(x) = A \x\~ 7 . These were determined by fitting the mean of the densities of 6 separate runs, over the range 
10~ 3 < x < 10 _1 , and bootstraping 10000 synthetic data sets to obtain the error bars. Note that there is some significant covariance between A and y. Column 
4 shows the reduced x^ ed - Column 5 shows the time average and variance of y for all time steps between 640 < t < 700 (see also Fig. pi. Columns 6 and 7 
show the inferred value of y obtained by fitting the cumulative mass (estimated from the data in two ways) over the same range in x, and Column 8 is inferred 
from a fit to the J(ju) relation in equation fl3) . All fits were performed at the end of our simulation, at time 700, except the J(ju) fit, which has to be done at 
much earlier times because the relation disintegrates (see section |4~5) . 



roughly unity if the data are consistent with the fitted density profile 
and the errors are Gaussian. However, we expect to be less than 
1 because the errors in the density bins show significant negative 
kurtosis. The apparent trend in y and A with N is illusory; it is 
only the case for this particular snapshot. Although the smallness 
of the bootstrap errors indicate that the fit is very stable at a given 
snapshot in time, the value of the best fit y varies from one snapshot 
to the next. 

Fig.[3]plots the best fit value of y for several snapshots with er- 
rors from the bootstrapped data, and demonstrates the extent of the 
time variability of the fit. Simulations with more particles display 
smaller absolute rms variability with time, but have a much larger 
disparity between the rms variability in time and the bootstrapped 
errors in a particular time step. This is because we have used a den- 
sity estimator with approximately 20 particles per bin for all reso- 
lutions in this study. Therefore the simulations with more particles 
have more data points available for the bootstrap. For all resolu- 
tions, however, the bin locations shift from one time step to the 
next, keeping the number of particles per bin roughly equal to 20. 
This procedure generates variability in the density estimator from 
one snapshot to the next. The time average and variance of y over 
the interval 640 < t < 700 are shown in column 5 of Table [T] 

As discussed in section [2] a self-similar solution exhibits its 



power-law behavior in several different functional relations. In 
principle, the exponent y can be extracted from measurements of 
any one of these relations. For example, the cumulative mass in- 
side a position x is m(x) oc \x\ l ~ 7 . Since the mass is a cumulative 
variable, mass data points are correlated; however we have opted 
to give all points equal weight in the bootstrap procedure. We have 
computed the cumulative mass in two ways, first by numerically 
integrating the mean density (abscissae at the bin locations), and 
second by avoiding the binning entirely and tallying the total mass 
at the particle locations for one of the runs (abscissae at the particle 
locations). The results of fitting the cumulative mass can be found 
in columns 6 and 7 of Table [T] The exponents are not consistent 
with the fits to the density, nor with one another. This discrepancy 
is not understood, though it is easy to identify systematics which 
may affect the fit. Being an integrated quantity, the cumulative mass 
m{x) is sensitive to all regions of the data interior to x; in particular, 
it is not possible to excise the innermost part of the object as we 
have done with the density. Including more of the central region in 
the density fits will cause the power-law exponent to be systemat- 
ically shallower than when it is excised, consistent with the mass 
result. The binning also likely affects the fit differently for the mass 
and the density. 

The most interesting of the indirect scaling relations is that 



© 0000 RAS, MNRAS 000, 000-000 



Gravitational Collapse in One Dimension 7 



0.001 



> 0.000 



-0.001 




Figure 4. Four variations in the initial velocity perturbation. The den- 
sity profiles arising from these initial states are shown in Fig. [5] and 
an intermediate snapshot of the phase space configuration is shown in 
Fig. [6] The functional forms are: [IC 0] i/ f (jcf) = -0.001 sin(x f ), [IC 1] 
Viixd = -0.0005 sin(xi), [IC 2] Vi (xi) = -0.001 sin(1.5 * Xi ), [IC 3] 
Vi (xi) = -0.001 sin 3 (jc f ) and [IC 4] i/ f (jc/) = +0.0001 sin(jc f ). 



expressed in equation p3) . The action of a particle orbit in the sta- 
tionary mean-field potential can be computed directly by integra- 
tion over the orbit. The values for exponents obtained by fitting this 
relation are shown in column 8 of Table Q] and show a remarkable 
stability as a function of N, compared to the other methods. In Fig. 
[2] the density profile corresponding to this value of the exponent is 
indicated by the dashed line. The discrepancy between the direct fit 
to the density and the indirect fit to the J(jj) relation is very clear 
in that figure. The details of the measurement and fit of the J(ji) 
relation are discussed further in section 03] 



4.2.3 Sensitivity to initial conditions 

To close this section on simulations of cold collapse, we discuss 
the sensitivity of the exponent y to changes in the initial conditions 
(ICs). Fig. [4] shows several of the initial conditions that we have 
tested. These are all cold and homogeneous in density, but with a 
variety of initial velocity perturbations. IC is the fiducial case, 
with particle velocities given by equation (T6] > with V = 0.001. IC 
1 tests the sensitivity to the velocity amplitude; IC 2 investigates 
the impact of a slow group of outer particles whose infall velocity 
lags that of intermediate particles; IC 3 has two additional inflec- 
tion points away from the origin; IC 4 examines initially expanding 
rather than initially contracting particles. 

Fig. [5] shows the final configuration (t = 700) of each IC. The 
top panel shows the average density of six simulations with the 
same number of particles as for the fiducial simulation, for each 
initial condition. The normalizations of all ICs except IC have 
been artificially shifted for visual clarity. Solid black curves with 
power-law behavior x~ 1/2 are also shown for reference. The bottom 
panel shows the same average densities, but with an x~ 1/2 trend 
scaled out. The normalizations are not shifted in this panel. We see 
that the resulting density profiles for these different initial condi- 
tions are consistent with an attractor behavior p(x) oc \x\ 7crit with 
Tcrit - 0.47. Power laws with the best fit y (in the legend) are plot- 
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Figure 5. Top: The density profiles at t = 700 for the initial conditions 
plotted in Fig. [4] These curves are averages over 6 simulations of ~ 10000 
particles, with slightly different particle numbers. The curves have been 
shifted vertically for visual clarity. The two solid black lines show x~ 1 ^ 2 
power laws for comparison. Bottom: The same densities scaled by x 1 ^ 2 , 
and not artificially shifted. This panel demonstrates more clearly that the 
attractor solution is somewhat shallower that p oc x -1 / 2 . The short verti- 
cal lines indicate the range of data that were used to fit the exponent and 
amplitude of the density. The lower curves show the best fit power law for 
each initial condition, also indicated in the legend. This helps quantify the 
extent to which the density approaches an attractor. We caution the reader 
that there is some covariance between the slope and the amplitude in the fit. 



ted in the lower part of the bottom panel, to help quantify the ex- 
tent to which the density approaches an attractor. The vertical black 
hashes show the range of data that was used in the fit. There is some 
covariance between the slope and the amplitude. 

Even though the final states of these simulations are similar, 
there are large qualitative differences in the formation history. Fig. 
[6] shows an intermediate snapshot of the phase-space configuration 
for these four simulations; ICs 1 and 2 undergo a monolithic col- 
lapse, while ICs 3 and 4 each form two distinct halos that later 
merge to form the final product. These can be compared with the 
bottom left panel of Fig.[T] which shows IC at the same time. The 
colors in Fig. [6] encode the initial positions of the particles; blue 
particles were initially located close to x = and red were at the 
largest distances. Note that the inner particles remain in the interior 
of the final halo only for cases that undergo monolithic collapse. 



4.3 Warm initial conditions 

The self-similar final states discussed above arise from cold ini- 
tial conditions. At the end of Section [2] we described "coldness" 
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Figure 6. An intermediate snapshot at time t = 40 in the evolution of ICs 
1-4. ICs 1 and 2 undergo monolithic collapse while ICs 3 and 4 first form 
two distinct objects that later merge to form the final virialized object. The 
color coding indicates the initial positions of particles. Blue particles started 
at small distance and red started at large distance. 
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Figure 7. A short wavelength component in the velocity perturbation 
makes the final density profile depart from a power law behavior. The cap- 
tion indicates the number of wavelengths over the interval (-n/2, n/2) in the 
initial conditions. As the wavelength approaches the Nyquist wavelength, 
equal to twice the initial inter-particle spacing, the simulation behaves more 
and more like a warm simulation (magenta curve) in which the initial veloc- 
ities have been perturbed with 5% random noise. These curves are averages 
over 6 simulations, and have been further averaged over 20 time units. The 
time averaging is used to smooth the impact of substructure on the density 
profile. The short wavelength component generates long lived substructure 
(cf. Fig. [8} that does not disappear over the timescales that we ran the sim- 
ulations. 



as localization of the particle density to a single smooth curve in 
the (x, z/)-plane. In this section we describe simulations whose fi- 
nal states are quite different from the self-similar states described 
above. We assume in this discussion that a state initialized with ve- 
locity perturbations from a stochastic process (noise) satisfies any 
conceivable definition of "warm". Our interest is in the approach to 
warmness as structure is added to the function describing the cold 
initial particle locus. We shall find that even a small amount of ad- 
ditional structure in the initial state, that appears consistent with the 
notion of "coldness," destroys the self- similarity of the final state. 

In a particle representation, the degree of smoothness of the 
particle locus must be compared to the particle separation. In the 
cold scenarios above, initial velocity perturbations were chosen to 
have wavelengths much larger than the initial particle separation. 
As the wavelength of the initial velocity perturbation decreases, it 
approaches the Nyquist limit for particle sampling, becoming more 
like a velocity noise in character. In this limit, we expect the subse- 
quent evolution to resemble that of warm initial conditions, where 
the velocity perturbations are generated by independent stochastic 
processes for each particle. As a consequence, we expect that the 
resulting equilibrium halo will develop a core; collapse from a state 
of nonzero velocity dispersion should develop a core because the 
maximum phase-space density cannot increase ( Tremai ne & Gunn| 
[1979). 

Fig. [7] shows the density profiles that result from adding an 
extra component to the velocity perturbation of equation (T6| . The 
additional component has a lower amplitude than the original per- 
turbation (0.1 V X but a higher wavenumber modulation. We have 
tested components containing 5, 50, 500 and 3000 wavelengths 
over the domain of the data; for comparison, the Nyquist rate (de- 
termined by the initial inter-particle spacing) would correspond to 
N/2 = 5000 wavelengths over the domain of the initial particle dis- 
tribution (-7r/2 < x < 7i /2). The bottom curve (magenta) shows 



a warm simulation in which an RMS 5% random fluctuation has 
been added to the initial velocities. These curves are averaged over 
six simulations and 20 time units (in the range t = 680-700) for 
reasons explained below. As expected, the density profile becomes 
shallower near the center as the spatial frequency of the initial ve- 
locity perturbation increases. These features indicate that it is dif- 
ficult to achieve coldness in particle realizations of the continuum. 
It is remarkable that a well-defined constant-density core appears 
even in a simulation for which the spatial frequency is only 10% of 
the Nyquist frequency. 

Less obviously, a high frequency perturbation in the initial ve- 
locities leads to the formation of substructures. This is illustrated 
in Fig. [8] which shows three successive times from the simulation 
plotted in blue of Fig. [7] in which the initial velocity perturbation 
contains five wavelengths in the domain of the initial particle dis- 
tribution. Many individual sub-halos form in locations where the 
flow is initially converging, before merging to form the final prod- 
uct. The sub-halos are very long lived; some of them never dissolve 
over the simulation times we have tested. Thus we can only mea- 
sure a density profile in a time averaged sense — this is why 20 time 
units have been averaged to produce the curves in Fig. [7] 

4.4 Cold initial conditions with inhomogeneous density 

All simulations discussed thus far were initialized with constant 
density. In this section we report on simulations initialized with an 
inhomogeneous density. In particular, we choose an initial power- 
law density profile p m [ t (x) oc \x\~ 7i . We will demonstrate an ap- 
parent bifurcation in the form of the final state as a function of 
the initial exponent y,-. The simulations described here all had 
N = 10006, 10012, 10018, 10024, 10030 and 10036 particles, 
with particle spacing chosen so that the resulting density is a power 
law with exponent y-. 
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Figure 8. Three snapshots from the simulation leading to the top (blue) density profile in Fig. [7] The small velocity perturbation results in the formation of 
several distinct objects that later merge to form the final product. The average central density at late times is less than that in simulations without small scale 
power. 
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Figure 9. Top: The final density profiles for initial conditions that follow 
a power law p oc \x\~ 7i . The final power law y is also noted in the caption. 
The curves have been artificially shifted for visual clarity. The solid black 
lines show |x| -1 / 2 power law for reference. Bottom: The same densities 
scaled by x 1/2 , and not artificially shifted. Initial power-law profiles that 
are shallower than p(x) oc |x| -1 / 2 evolve toward an approximate attractor, 
but initial profiles that are steeper, such as the red curve, retain their initial 
power law. 



Fig. [9] summarizes the results of a set of simulations with ini- 
tial exponents y f = 0.8,0.5,0.3, -0.5. The curve labeled y- = 0.0 
shows our fiducial simulation for comparison. All initial particle 
velocities were set to v = in these simulations (except the fiducial 
case). We find that initial conditions with slope y t < y crit collapse 



to form final states with roughly the same \x\~ 7cTit density profile as 
homogeneous initial states, shown in dark blue in Fig. [9] However, 
initial conditions with y- > y cr i t , such as the red curve in Fig. |9| all 
preserve their density profiles. The bottom panel of Fig. [9] shows 
that the attractor is not as tight as in the homogeneous case; the 
curves do not overlay one another as tightly as in Fig. [5] and there 
is more scatter in the measured power law index. Also the central 
portions of the profile appear to be flatter than the behavior at larger 
|jc|. For this reason the fit values (displayed in the legend) are some- 
what more sensitive than fits in Fig. [5] to the range of data used in 
the fit (indicated in the bottom panel of Fig. [9] with short vertical 
black lines). 

Fig. 1 10| shows the exponent of the final density profile, y, ver- 
sus that of the initial condition for a larger set of values of y-. 
Each data point is computed using 6 simulations. For these simula- 
tions, we have perturbed the particle velocities according to equa- 
tion (T6] >, because this significantly decreases the run time. For the 
cases with y = 0.8, 0.5, 0.3, -0.5, we have verified that the veloc- 
ity perturbation does not significantly impact the measured slopes. 
Solid lines for y = y- and y = y crit are also plotted for reference. 

Although the bootstrap method significantly reduces the sen- 
sitivity to the range of data used for the fits in Fig. [To] it does not 
eliminate it entirely. If either the inner or outer boundary of the 
domain is reduced by an order of magnitude or more, the best fit 
y becomes systematically shallower. The fitted exponents in Fig. 
[To] are power-law fits to the portion of the density in the range 
10~ 3 < |x| < 10 _1 . The figure suggests that a system initialized in 
a more collapsed state than the attractor solution remains close to 
its initial power law. However, a system initialized with a shallower 
density profile than the attractor significantly steepens to evolve to- 
ward the attractor. The dip of very shallow slopes for < y < 0.2 
has its origin in two features of these final states. First, at t = 200 
the simulations with < y < 0.2 display noticeable gaps in phase 
space structure in the range 10~ 3 < \x\ < 10" 1 . These gaps, and 
similar gaps in the homogeneous simulations, are discussed below. 
These gaps affect the fit for the exponent. Second, simulations with 
< y < 0.2 achieve much lower peak phase-space densities than 
the homogeneous case or cases with y > 0.2. We postulate that 
either the gaps or other discreteness effects interfere with the effi- 
cient transfer of energy from the inner particles to the outer ones. 
Since the inner threshold is fixed, we are not discarding the very in- 
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Figure 10. The final exponent y of the density profile p(x) oc \x\~ r at time 
t = 200, as a function of the initial exponent y, in the initial inhomogeneous 
distribution of particles Pinit(^) 00 \A~ 7i • Also shown in green are fits using 
the integrated mass interior to \x\. Initial slopes steeper than the attractor 
solution j{ > y cr it tend to steepen only slightly. In contrast, initial slopes 
much shallower than the attractor yi < y cr jt steepen considerably. 

nermost portion of the data in these simulations, and the fits result 
in systematically shallower values of y, as discussed above. 

4.5 Time evolution of rank, action, energy and density 

Fig. [TT] shows the evolution of the particle energies and actions as 
a function of time. The different colors mark quintiles in the initial 
particle position, as described in Section]!] The evolution of the en- 
ergies for the fiducial simulation, in the top left corner, shows rapid 
oscillations for t < 50. The phase space plots of the system config- 
uration can be seen in Fig. [13] at the times indicated by the vertical 
magenta lines. Before time 22, the system is undergoing rod-like 
rotation in phase space, with only one stream per location x. The 
oscillations in the energy are dominated by changes in the potential, 
as the system oscillates between maximum rarefaction and maxi- 
mum compression. At approximately t = 22, the system transitions 
from one to three streams at small x. After this time the central part 
of the particle locus winds up quickly, with more than 15 streams 
in the inner region by time t = 25, but the oscillation pattern in 
energy continues in phase with the outer particles, suggesting that 
the potential in the center is dominated by the configuration of the 
particles in the outskirts at these times. The range of particle ener- 
gies spreads with time, indicating that energy is being transferred 
from inner to outer particles. The colors are not mixed, demonstrat- 
ing that the particles approximately maintain their rank ordering in 
energy; violent relaxation redistributes energy among the particles 
but does not scramble the energy ordering. For t > 50 the ener- 
gies of the particles are nearly constant. For example, an extended 
run shows that the energies of the (most, least) energetic particles 
change by (0.5%, 1.0%) between 150 and 700 time units. This sim- 
ulation was run for several values of N, as indicated in Table [T] The 
energy oscillations from violent relaxation occur in phase for all 
of these runs, and the final energies of the (innermost, outermost) 
particles of the two highest resolution runs match to within (0.7%, 
0.6%) in all the simulations. 

The time evolution of the particle actions in the fiducial simu- 



lation is shown in the bottom left hand panel of Fig.[TT] The actions 
are computed assuming the potential is fixed at its instantaneous 
functional form, even though for times less than t = 50 the po- 
tential is evolving significantly. The fluctuations in low / at late 
times are almost certainly caused by discreteness noise. The ac- 
tions display a number of interesting features. First, they stabilize 
extremely abruptly, well before the energy oscillations damp out 
and the energies of the particles cease to evolve. Second, the transi- 
tion to constant action for the inner particles happens just after the 
transition from rod-like rotation to multiple stream flow. This can 
be seen in the fiducial simulation by comparing to the top row of 
Fig. [13] which shows that the transition from rod-like rotation to 
multi- stream flow happens at approximately time t = 22. 

Comparing the energy and action plots we learn that the re- 
laxation of the system appears to be a two stage process. The first 
stage, in which both the energies and the actions evolve rapidly, 
happens primarily before the transition from one to three streams. 
The second stage, in which the actions rapidly settle and the en- 
ergies evolve adiabatically, happens as the system winds up, pri- 
marily between the one-to-three stream transition, and the ultimate 
stabilization of the outer gravitational potential. During this sec- 
ond stage, the amount of energy transferred between particles, and 
the distances over which it can be transferred both tend to zero as 
the number of streams increases and the inter-caustic distance de- 
creases. 

The panels on the right of Fig. show the evolution of the 
particle energies and actions for one of the simulations initialized 
with a steep power law density profile. These plots contrast starkly 
with the evolution of the fiducial simulation. We see that the ac- 
tions (bottom right panel) are immediately stable, and evolve very 
little from their initial values. This is the case because the one-three 
stream transition occurs well before time t = 1, and the system does 
not undergo several rod-like rotations before the transition to mul- 
tiple streams. Therefore, steep systems like this one relax mainly 
through the second adiabatic stage of energy redistribution. The re- 
distribution of particle energies (top right panel) is much less than 
in the fiducial case. Since the particles are essentially in fixed sta- 
ble orbits from time t = 1 onward, the change in the scaling of the 
density profile is small, as can be seen in Fig.fTo] 

A curious feature in the top left panel of Fig.[TT]is the devel- 
opment of narrow gaps in the energy distribution, seen after about 
60 time units, in the yellow and green curves. The gaps correspond 
to specific phase space pockets, empty of particles, which open at 
times significantly after the end of violent relaxation. An example 
can be seen in the top panel of Fig. [12] The emergence of these 
gaps scrambles the energy rank-ordering of the particles that had 
been preserved in the earlier stages of the evolution (bottom panel). 

These gaps are a robust result of our simulations. We observe 
them in simulations with all values of N, with and without impos- 
ing explicit left/right symmetry, and gaps appear in both the heap 
and the search-tree versions of the code run on identical initial con- 
ditions. They also occur in the simulations initialized with a power 
law density profile, though they manifest at different times and are 
much less pronounced when the initial density profile is steep. 

Fig. [13] shows the state of the fiducial system at 8 specific 
times, which are indicated in the left panels of Fig. [TT] with ver- 
tical lines. The phase- space distribution (first row) is shown with 
color indicating quintiles in initial position. The evolution of the 
density profile (second row) shows the steady growth in the num- 
ber of caustics. At early times the density is roughly constant in 
the intervals between caustics. As the evolution proceeds, the caus- 
tics become more closely spaced and the usual (x max - x)~ l/2 den- 
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Figure 11. A comparison between the fiducial simulation, initialized with a homogeneous density, and a simulation in which the initial configuration of 
particles follows a steep power law density profile. Top row: The evolution of the energy of every 40th particle in the fiducial simulation, and the simulation 
with a steep initial power law configuration with y; = 0.8. The period of violent relaxation in the fiducial simulation can be seen between times and -50. 
After this time the gravitational potential is approximately stationary and the particle energies hardly evolve. The vertical lines indicate times shown in Fig. 

t Bottom row: The actions for the fiducial and power law simulations, derived by assuming a stationary potential and performing the integral of equation 
. The fluctuations for low J at late times are almost certainly caused by discreteness noise. 



sity distribution associated with a fold singularity at x max dominates 
the behavior between caustics. We also notice that the caustics are 
roughly equally spaced in logarithmic distance. 

Also shown in Fig. [13] are plots of energy versus position for 
the particles (third row). The colors in this row again indicate quin- 
tiles in initial particle position. As seen in these plots, violent re- 
laxation mixes the final positions relative to the initial ones (the 
colors are mixed when projected onto the horizontal axis) but does 
not mix the final energies relative to the initial ones (the colors are 
not mixed when projected onto the vertical axis). The envelope of 
the energy-position plot gives the maximum orbital excursion as a 
function of energy; the locus of this envelope is the potential as a 



function of distance, and therefore displays the power-law behavior 
| x |2-rcrit p or comparison, this form of the potential is shown with a 
green line. Indeed, if the turnaround is known for all the particles, 
this behavior can be used as another method to infer y. We exam- 
ined this for the fiducial simulations and found that the power law 
is slightly shallower than when fitting the cumulative mass, though 
not as shallow as when fitting J (pi). Since the potential is a double 
integral of the density, it displays correlations and dependency on 
the global structure of the state, similar to the effects discussed for 
the cumulative mass scaling in section [42] 

The evolution of the J(ja) relation is displayed in Fig. [T3| 
(fourth row) and shows a very clean power law developing even be- 
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Figure 12. Top: The gaps developing in the energy distribution in Fig. 
[TT] correspond to empty pockets in phase space. This panel shows gaps in 
the yellow and green marked particles, at energies E - O m i n « 0.6 and 
E - Omin ~ 0-9. Bottom: The energy gaps cause the rank ordering of the 
particles to be scrambled, in certain bands. The insets show the rank order- 
ing at two earlier times, also studied in Fig.[T3] 



fore the end of violent relaxation (t = 25). Here we have scaled by 
IT 2 so that the initial relation is flat. The solid black line shows the 
scaling from equation fl3) for the value y crit = 0.47, demonstrat- 
ing that the asymptotic scaling of this relation predicts a different 
value of y than the other scaling relations. The power-law behavior 
develops much earlier than the time at which power-law behavior 
in the density becomes apparent. The panel at t = 22 shows that 
this power-law behavior is first seen in the outer parts of the system 
and last established in the inner parts, which is perhaps somewhat 
counterintuitive. This dynamical behavior is related to the fact that 
the phase-line remains close to a straight line longest at x = 0. The 
clean power law in JQi) deteriorates at later times, apparently due 
to the loss of coherence in the center of the system, consistent with 
the loss of correlation between fi and particle energy. 

As already indicated in Table [T] the value of the exponent y 
extracted from the measured J(jj) relation is systematically lower 
than that extracted from fits to the density. However, these expo- 



nents cannot be compared at the same times, because the J(p) rela- 
tion deteriorates before the density is sufficiently settled. All fits to 
the J(ja) relation referred to in this paper have been done at times 
around t = 30 to allow sufficient dynamic range for an accurate fit, 
10~ 2 < \i < 10 _1 . However, the observed slope of J{ji) does not ap- 
pear to evolve at all past this time, and it is easy to conjecture that a 
simulation with sufficient resolution in the center to delay the dete- 
rioration would show no subsequent change or drift in the exponent 
determined from J(jS). We conjecture that the reason the exponent 
obtained from J(ja) does not agree with fits to the density is because 
the monotonicity of the relation is destroyed by discreteness effects 
(see bottom panel of Fig. [TTJ well before the system is sufficiently 
wound to render the expression in equation fT3] ) valid. 

The bottom row in Fig. [13] shows the evolving distributions 
of orbital actions (per unit we expect a coarse-grained power- 
law behavior as expressed in equation (T2] >. This is shown in the 
plot with a solid black line for y crit = 0.47. The observed distri- 
bution shows clear systematic deviations from the expected power 
law at intermediate times and becomes sufficiently disordered at 
later times that it is impossible to extract a meaningful value of y. 
Although we might expect certain subtleties in the convergence of 
the distribution in /, it is somewhat surprising that this distribu- 
tion can display such complex structure, while the observed J{p) 
relation is quite stable and clean. 

We do not show similar snapshots for the simulations initial- 
ized with a power-law density profile; however, we now describe 
the results qualitatively. The violent relaxation and settling to an 
equilibrium gravitational potential proceed much faster in simula- 
tions with a steep initial profile compared to those with a shallow 
initial profile. The development of gaps of the type shown in Fig. [T2] 
is strongly suppressed in the steep cases, compared to the shallow 
cases. Finally, there is a difference in the evolution of the distribu- 
tion of actions: the shallow cases eventually reach an f(J) distribu- 
tion that follows equation |T2| with respect to the final slopes of the 
density profiles. However, the steep cases do not. The f(J) distribu- 
tion for the steep cases begins very close to the f(J) that would be 
expected for a p oc x~ 7ciit scaling, and since the actions retain their 
initial values (shown in Fig. [TTJ, this distribution does not evolve at 
all. 



5 CONCLUSIONS 

We have performed a suite of simulations of one dimensional grav- 
itational collapse, equivalent to N infinite parallel mass sheets col- 
lapsing along one axis in three dimensions. We find that a self- 
similar attractor solution exists for "cold" initial conditions with 
near-homogeneous density; the density in the final state attractor 
follows a power law, p oc |x| _rcrit , over at least three decades of 
distance, with y crit 0.47. The measured value of the exponent is 
shallower than the value 1/2 conjectured by Binney (2004). This 
self- similar state arises after a period of violent relaxation that re- 
distributes energy while preserving the energy ranking of the parti- 
cles. The measurements of the exponent for the power-law density 
profile are done at times significantly after the end of violent re- 
laxation (t = 700 versus t = 60), to ensure that the coarse-grained 
density has achieved a stationary form. Measurements of the ex- 
ponent obtained from the J(ji) relation can (and must) be done at 
significantly earlier times. These show that an attractor solution is 
in place at very early times. However, the exponent obtained from 
J(jj) is systematically lower than that obtained from the density at 
later times. The two determinations cannot be compared because at 
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Figure 13. Eight snapshots of the fiducial simulation are shown. The time is given at the top of each column; each time is indicated in Fig.^^by a vertical 
line. The top row shows the phase-space distribution of the particles. The colors indicate 5 quintiles in initial particle position, blue started near x = and 
red started far away. The second row shows the density profile. The diagonal black line is x~ 1 ^ 2 . The third row shows a scatter plot of the energy and current 
position, the colors again indicating quintiles in initial position. For density p oc |x| -7crit , the potential varies as |x| 2-7crit and this is shown with a solid green 
line using y cr i t = 0.47. The fourth row shows the measured J(p) relation (equation [T3). The bottom row shows the distribution of the number of particles per 
unit /. 
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later times discreteness effects destroy the phase coherence of the 
simulations, and with it the monotonicity of J(ji). If the initial con- 
ditions are warmed by applying smooth velocity perturbations with 
spatial frequencies approaching the Nyquist limit, a cross-over is 
reached where self- similarity of the final state is destroyed and is 
replaced by a cored solution. Cold initial conditions with a steep 
initial power-law density of exponent y t > y crit do not reach the at- 
tractor either; they instead collapse to form a final equilibrium state 
with 7 ~ j[. 

Our simulations provide some insight into the process of vi- 
olent relaxation and provide hints that may be useful for the de- 
velopment of an analytical understanding of the properties of the 
attractor solution. For example, we find that many of the scaling 
laws are in place very early, even during the violent relaxation pro- 
cess. The evolution of some of the scaling relations suggests that 
the violent relaxation process finishes first in the outer parts of the 
structure, while the inner particles take longer to redistribute their 
energies. Such insights may then also be extended to the 3D case, 
eventually shedding light on the origin of cold dark matter halo 
density profiles. 
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